use "HSSC replication data\Gallup data\cleaned_worldpoll_microdata.dta", clear

**Google mobility cumulative-to-date
merge m:1 countrycode month year using  "HSSC replication data\Public data\Google_Mobility_Data_March31_2021_by_Month_Year_cum_to_date.dta"
tab country if _merge==1
drop if _merge==2
drop _merge

**Merge cumulative-to-date stringency data
merge m:1 countrycode month year using "HSSC replication data\Public data\Oxford_Policy_Data_Month_Cumulative_to_Date.dta"
drop if _merge==2
drop _merge

*IHME cumulative-to-date
merge m:1 countrycode month year using "HSSC replication data\Public data\IHME_excess_deaths_monthly_cum_to_date.dta"
drop if _merge==2
drop _merge

label var employed "Worked for pay duing past 7 days"
label var RTstringencyindex "Stringency index--cumulative-to-date (standardized)"
label var  RTeconomicsupportindex "Economic support index--cumulative-to-date (standardized)"
label var  RTIHME_reported_deaths "Officially reported COVID-19 deaths per capita"
label var RTIHME_est_deaths  "Estimated actual COVID-19 deaths per capita"
label var RTretail_7day "Visits to restaurants and related establishments"

label var incQ1 "Bottom income quintile (relative to top-fifth)"
label var incQ2 "Second income quintile (relative to top-fifth)"
label var incQ3 "Third income quintile (relative to top-fifth)"
label var incQ4 "Fourth income quintile (relative to top-fifth)"
label var city "City (relative to rural)"
label var suburb "Suburb (relative to rural)"
label var town "Town (relative to rural)"
label var elementary "Elementary education (relative to tertiary)"
label var secondary "Secondary education (relative to tertiary)"

label var harm_index"Harm index"
label var lost_job "Lost job"
label var lost_income "Lost income"
label var affected_alot "Affected a lot"
label var temp_layoff "Temporarily laid off"
label var lost_hours "Lost hours"

global restrictions "stringencyindex c1_schoolclosing c2_workplaceclosing c3_cancelpublicevents c4_restrictionsongatherings c5_closepublictransport c6_stayathomerequirements c7_restrictionsoninternal c8_internationaltravel h1_publicinformation h2_testingpolicy h3_contacttracing h4_emergency h5_investmentinvaccines h6_facialcoverings h7_vaccinationpolicy h8_protectionofelderlypeople"
global controls "out_labor city suburb town foreign elementary secondary ageg1 ageg2 ageg3 ageg4 ageg5 ageg6   married have_children mother female  incQ1 incQ2 incQ3 incQ4 "

********
****SUMMARY STATISTICS
****
sum RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day [aw=WGT], sep(0)

sum RTstringencyindex RTc1_schoolclosing RTc2_workplaceclosing RTc3_cancelpublicevents RTc4_restrictionsongatherings RTc5_closepublictransport RTc6_stayathomerequirements RTc7_restrictionsoninternal RTc8_internationaltravel RTh1_publicinformation RTh2_testingpolicy RTh3_contacttracing RTh4_emergency RTh5_investmentinvaccines RTh6_facialcoverings RTh7_vaccinationpolicy RTh8_protectionofelderlypeople RTeconomicsupportindex ///
RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day , sep(0)

sum RTstringencyindex RTc1_schoolclosing RTc2_workplaceclosing RTc3_cancelpublicevents RTc4_restrictionsongatherings RTc5_closepublictransport RTc6_stayathomerequirements RTc7_restrictionsoninternal RTc8_internationaltravel RTh1_publicinformation RTh2_testingpolicy RTh3_contacttracing RTh4_emergency RTh5_investmentinvaccines RTh6_facialcoverings RTh7_vaccinationpolicy RTh8_protectionofelderlypeople RTeconomicsupportindex ///
RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day if month>=3 & month<=5 & year==2020, sep(0)

sum RTstringencyindex RTc1_schoolclosing RTc2_workplaceclosing RTc3_cancelpublicevents RTc4_restrictionsongatherings RTc5_closepublictransport RTc6_stayathomerequirements RTc7_restrictionsoninternal RTc8_internationaltravel RTh1_publicinformation RTh2_testingpolicy RTh3_contacttracing RTh4_emergency RTh5_investmentinvaccines RTh6_facialcoverings RTh7_vaccinationpolicy RTh8_protectionofelderlypeople RTeconomicsupportindex ///
RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day if month>=6 & month<=9 & year==2020, sep(0)

sum RTstringencyindex RTc1_schoolclosing RTc2_workplaceclosing RTc3_cancelpublicevents RTc4_restrictionsongatherings RTc5_closepublictransport RTc6_stayathomerequirements RTc7_restrictionsoninternal RTc8_internationaltravel RTh1_publicinformation RTh2_testingpolicy RTh3_contacttracing RTh4_emergency RTh5_investmentinvaccines RTh6_facialcoverings RTh7_vaccinationpolicy RTh8_protectionofelderlypeople RTeconomicsupportindex ///
RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day if month>=10 & month<=12 & year==2020, sep(0)

sum RTstringencyindex RTc1_schoolclosing RTc2_workplaceclosing RTc3_cancelpublicevents RTc4_restrictionsongatherings RTc5_closepublictransport RTc6_stayathomerequirements RTc7_restrictionsoninternal RTc8_internationaltravel RTh1_publicinformation RTh2_testingpolicy RTh3_contacttracing RTh4_emergency RTh5_investmentinvaccines RTh6_facialcoverings RTh7_vaccinationpolicy RTh8_protectionofelderlypeople RTeconomicsupportindex ///
RTIHME_reported_deaths RTIHME_est_deaths  RTretail_7day if month>=1 & month<=3 & year==2021, sep(0)


********
**Main models
***********
eststo clear

foreach Y in harm_index  lost_job lost_income affected_alot temp_layoff lost_hours  {
eststo: quietly mixed `Y'  RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths RTretail_7day $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths   $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
eststo: quietly mixed `Y'  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths  RTretail_7day $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
}

foreach Y in harm_index  lost_job lost_income affected_alot temp_layoff lost_hours  {
eststo: quietly mixed `Y'  Quin_stringency1 Quin_stringency2 Quin_stringency3 Quin_stringency4 RTeconomicsupportindex RTIHME_est_deaths  RTretail_7day $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
}

esttab using "HSSC replication data\Findings\multi_level_model_st_errors.csv", ///
stats(log_likelihood N, label("Log likelihood" "Sample Size")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01)   nogaps replace


**ALT HARM***
foreach x in harm_index harm_pca covid_effect {
egen Z`x'=std(`x')
}
eststo: mixed Zharm_index  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths   $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
eststo: mixed Zharm_pca  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths   $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
eststo: mixed Zcovid_effect  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths   $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
eststo: mixed Zharm_index  RTstringencyindex RTeconomicsupportindex RTIHME_est_deaths   $controls if SDtemp_layoff!=. & SDlost_job!=. & SDlost_hours!=. & SDlost_income!=. & SDaffected_alot!=. [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
esttab using "HSSC replication data\Findings\Sup_multi_level_model_alternative_harm_index.csv", ///
stats(log_likelihood N, label("Log likelihood" "Sample Size")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01)   nogaps replace


***********
**Policy decomposition
***********
eststo clear

foreach X in  stringencyindex c1_schoolclosing c2_workplaceclosing c3_cancelpublicevents c4_restrictionsongatherings c5_closepublictransport c6_stayathomerequirements c7_restrictionsoninternal c8_internationaltravel h1_publicinformation h2_testingpolicy h3_contacttracing h4_emergency h5_investmentinvaccines h6_facialcoverings h7_vaccinationpolicy h8_protectionofelderlypeople {
eststo: quietly  mixed harm_index  RT`X' RTeconomicsupportindex RTIHME_reported_deaths  $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
gen B_`X'=_b[RT`X']
gen SE_`X'=_se[RT`X']
}

foreach X in  stringencyindex c1_schoolclosing c2_workplaceclosing c3_cancelpublicevents c4_restrictionsongatherings c5_closepublictransport c6_stayathomerequirements c7_restrictionsoninternal c8_internationaltravel h1_publicinformation h2_testingpolicy h3_contacttracing h4_emergency h5_investmentinvaccines h6_facialcoverings h7_vaccinationpolicy h8_protectionofelderlypeople {
eststo: quietly  mixed lost_job  RT`X' RTeconomicsupportindex RTIHME_reported_deaths  $controls [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  
estadd local log_likelihood `e(ll)', replace
gen B2_`X'=_b[RT`X']
gen SE2_`X'=_se[RT`X']
}

esttab using "HSSC replication data\Findings\multi_level_model_by_policy_st_error.csv", ///
stats(log_likelihood N, label("Log likelihood" "Sample Size")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01)   nogaps replace


***
**Interactions with income
******
global interact_controls "out_labor city suburb town foreign ageg1 ageg2 ageg3 ageg4 ageg5 ageg6   married have_children mother female"
eststo clear


foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
gen Stx_`x'=`x'*RTstringencyindex
}

eststo: quietly mixed harm_index  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4  RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Pharm`x' if `x'==1
}

eststo: quietly mixed lost_job  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4 RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Pjob`x' if `x'==1
}

eststo: quietly mixed affected_alot  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4 RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Paffect`x' if `x'==1
}

eststo: quietly mixed temp_layoff  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4 RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Ptemp`x' if `x'==1
}

eststo: quietly mixed lost_income  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4 RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Pinc`x' if `x'==1
}

eststo: quietly mixed lost_hours  Stx_incQ1 Stx_incQ2 Stx_incQ3 Stx_incQ4 incQ1 incQ2 incQ3 incQ4 RTstringencyindex RTeconomicsupportindex RTIHME_reported_deaths  $interact_controls  [pw=WGT] || WP5: || yrm:, mle nolog vce(cluster WP5)  

foreach x in incQ1 incQ2 incQ3 incQ4 incQ5  {
predict Phours`x' if `x'==1
}

esttab using "HSSC replication data\Findings\multi_level_model_interactions_with_income_st_error.csv", ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01)   nogaps replace

**Policy heterogeneity
keep in 1
keep B2_* SE2_* B_stringencyindex SE_stringencyindex B_c1_schoolclosing SE_c1_schoolclosing B_c2_workplaceclosing SE_c2_workplaceclosing B_c3_cancelpublicevents SE_c3_cancelpublicevents B_c4_restrictionsongatherings SE_c4_restrictionsongatherings B_c5_closepublictransport SE_c5_closepublictransport B_c6_stayathomerequirements SE_c6_stayathomerequirements B_c7_restrictionsoninternal SE_c7_restrictionsoninternal B_c8_internationaltravel SE_c8_internationaltravel B_h1_publicinformation SE_h1_publicinformation B_h2_testingpolicy SE_h2_testingpolicy B_h3_contacttracing SE_h3_contacttracing B_h4_emergency SE_h4_emergency B_h5_investmentinvaccines SE_h5_investmentinvaccines B_h6_facialcoverings SE_h6_facialcoverings B_h7_vaccinationpolicy SE_h7_vaccinationpolicy B_h8_protectionofelderlypeople SE_h8_protectionofelderlypeople
gen n=1
reshape long B_ SE_ B2_ SE2_, i(n) j(policy) string

replace policy="stringency index" if policy=="stringencyindex"
replace policy="school closing" if policy=="c1_schoolclosing"
replace policy="workplace closing" if policy=="c2_workplaceclosing"
replace policy="cancel public events" if policy=="c3_cancelpublicevents"
replace policy="restrictions on gatherings" if policy=="c4_restrictionsongatherings"
replace policy="close public transport" if policy=="c5_closepublictransport"
replace policy="stay-at-home-order" if policy=="c6_stayathomerequirements"
replace policy="internal travel" if policy=="c7_restrictionsoninternal"
replace policy="international travel" if policy=="c8_internationaltravel"
replace policy="public information" if policy=="h1_publicinformation"
replace policy="testing policy" if policy=="h2_testingpolicy"
replace policy="contact tracing" if policy=="h3_contacttracing"
replace policy="emergency declaration" if policy=="h4_emergency"
replace policy="vaccine investment" if policy=="h5_investmentinvaccines"
replace policy="mask order" if policy=="h6_facialcoverings"
replace policy="vaccination policy" if policy=="h7_vaccinationpolicy"
replace policy="protection of elderly" if policy=="h8_protectionofelderlypeople"
tab policy

*Drop investment items
drop if policy=="emergency declaration"
drop if policy=="vaccine investment"

gen effect=B_*100
gen up=100*(B_+(1.96*SE_))
gen low=100*(B_-(1.96*SE_))

gen effect2=B2_*100
gen up2=100*(B2_+(1.96*SE2_))
gen low2=100*(B2_-(1.96*SE2_))

cd "HSSC replication data\Findings\"

*Parsimonious--harm
graph dot (mean)  effect up low  , over(policy, label(labsize(small)) sort(1) descending) ///
marker(1, msize(large) msymbol(Oh)  mcolor(blue*.7)) ///
marker(2, msize(small) msymbol(Dh) mcolor(red*.7)) ///
marker(3, msize(small) msymbol(Dh) mcolor(red*.7)) ///
yscale(range(-15 40))  ylabel(-15(5)40) ///
saving(policy1) ///
pcycle(5) ///
legend(order(1 "Mean effect in std dev" 2 "95% c.i." )  region(style(none))  rows(1) span)  ///
ytitle("Estimated effect", margin(medsmall)) ///
title("A. Predicted effect on harm index by policy", span) ///
plotregion(color(white)) graphregion(color(white)) yline(0)  
graph export "HSSC replication data\Findings\Fig_Harm_ML_Regress_on_Alt_Stringency.png", replace

*job loss
graph dot (mean)  effect2 up2 low2  , over(policy, label(labsize(small)) sort(1) descending) ///
marker(1, msize(large) msymbol(Oh)  mcolor(blue*.7)) ///
marker(2, msize(small) msymbol(Dh) mcolor(red*.7)) ///
marker(3, msize(small) msymbol(Dh) mcolor(red*.7)) ///
yscale(range(-15 30))  ylabel(-15(5)30) ///
saving(policy2) ///
pcycle(5) ///
legend(order(1 "Mean effect in std dev" 2 "95% c.i." )  region(style(none))  rows(1) span)  ///
ytitle("Estimated effect", margin(medsmall)) ///
title("B. Predicted effect on job loss by policy", span) ///
plotregion(color(white)) graphregion(color(white)) yline(0)  
graph export "HSSC replication data\Findings\Fig_job_loss_ML_Regress_on_Alt_Stringency.png", replace

save "HSSC replication data\Findings\harm_coef_ml_model_by_policy.dta", replace

cd "`root_dir'\HSSC replication data\Findings\"

grc1leg policy1.gph policy2.gph    ///
, imargin(.5 .5 .5 .5) row(2) col(1) legendfrom(policy1.gph) iscale(*.85) plotregion(color(white)) graphregion(color(white))  
graph export "Fig_Harm_and_Job_Loss_ML_Regress_on_Alt_Stringency.png", replace

